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ABSTRACT 

We study the spatial distribution of dark matter halos in the Universe in terms of their number 
density contrast, related to the underlying dark matter fluctuation via a non-local and non- linear bias 
random field. The description of the matter dynamics is simplified by adopting the 'truncated' Zel'dovich 
approximation to obtain both analytical results and simulated maps. The halo number density field in our 
maps and its probability distribution reproduce with excellent accuracy those of halos in a high-resolution 
A-body simulation with the same initial conditions. Our non-linear and non-local bias prescription 
matches the iV-body halo distribution better than any Eulcrian linear and local bias. 

Subject headings: galaxies: statistics - large-scale structure of Universe 



1. INTRODUCTION 

The simplest description for biasing assumes that the 
fluctuations in the number density of luminous objects, 
Siumy and in the mass, <5 mass , are proportional, S\ um = 
b5 mas8 , where b is the so called linear bias factor. Re- 
cently, Catelan et al. (1998, CLMP), following the seminal 
papers by Cole & Kaiser (1988) and Mo & White (1996, 
MW), showed how the relation between dark halos, recip- 
ient of the luminous matter, and the underlying mass is 
to be cast in terms of a bias random field b, which de- 
pends in a non-local way on the mass density field. Halo 
biasing is a process which evolves in time, depends on the 
scales and the collapse times of the selected objects, but is 
additionally determined by the gravitational conditions of 
the environment. Most important, unlike previous models, 
CLMP treated halo biasing as a non-local process. 

In this Letter, we apply the CLMP bias model to ana- 
lyze the spatial halo distribution at several scales. Mass 
particles move according to the Zel'dovich (1970) approx- 
imation. We generally find excellent agreement between 
our theoretical predictions and the distribution of halos 
extracted from an A-body simulation with the same ini- 
tial conditions. In § 2 we present our bias model, in § 3 we 
test it against simulations; §4 contains our conclusions. 

2. THE DISTRIBUTION OF HALO FLUCTUATIONS 

2.1. The model 

Let us consider a population of halos, selected in La- 
grangian space through their mass M and formation red- 
shift Zf. At any comoving position x and observation 
redshift z < Zf, one can generate Eulerian maps of the 
halo number density fluctuation 5/i(x, z\M, Zf), given the 
mass density contrast 5{x, z) with Lagrangian resolution 
R and the corresponding Lagrangian halo density fluctu- 
ation field, 6fc(q\M, Zf), through the relation 



6 h (x,z\M,z f ) = [l + S%(q\M,z f )] [1 + <5(x,z)] - 1 , (1) 



valid in the single-stream regime (CLMP). The non- 
locality comes from the fact that the halo number density 
in x is determined by its initial value at the Lagrangian 
comoving position q. Using a local version of the Press & 
Schechter (1974, PS) approach, we obtain 

8h(*l\M,Zf) = V^[tf- e (q)] 6 [tf - e (q)] 
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This result is obtained as follows: adopting the peak- 
background split, in the patch of fluid at q with La- 
grangian size R , one writes the mean number-density of 
halos in terms of the PS formula, with a local collapse 
threshold tf — e a modulated by the background field e , 
and finally removes the overall mean halo number density 
[eq.(3) below]. Here tf = 6 C /D(zf), with S c the critical 
threshold for collapse of a spherical perturbation and D(z) 
the linear growth factor of density fluctuations normalized 
to unity at z — [in the Einstein-de Sitter universe, S c ~ 
1.686 and D(z ) = (1 + z) -1 ]; e is the linear mass fluctua- 
tion extrapolated to z = and smoothed on R , with a 2 its 
variance. Finally, a 2 = (27T 2 )- 1 dk k 2 P{k)W{kR) 2 , 
with W(kR) the filter function and P(k) the primordial 
power spectrum, is the variance on scale M of the lin- 
ear density field e M and S 2 = a 2 — a 2 . Eq.(2) actu- 
ally generalizes eq.(42) of CLMP in that collapse on the 
background scale R is accounted for by the step func- 
tion <d[tf — e ]: halos of mass M oc R 3 cannot be present 
in a collapsed region of Lagrangian size R a > R. As 
stressed by CLMP (see also Porciani et al. 1998), this ap- 
proach defines catalogs of halos unaffected by the cloud- 
in-cloud problem (e.g. Bond et al. 1991) up to the scale 
R D . By expanding eq.(0) to first order in e one ob- 
tains, for do <C a M , 8%(q\M,Zf) ~ 6f(M, Zf) e Q (q), where 
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bi(M,Zf) = [tf/cr^ — l/tf] is the linear Lagrangian bias 
factor (MW). Note that in our approach the background 
scale R appears as a fitting parameter which can be used 
to optimize the performance of the model. As we will see 
in § 3.2, the optimal value of 0~ o generally depends on the 
chosen Eulerian resolution scale and on the masses of the 
considered halo population. 

When comparing with halos in numerical simulations, 
we will consider finite mass intervals, so we will have to 
replace 5^ in eq.(2) with its weighted average, where the 
weight is given by the comoving conditional mass-function 
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with pi, the mean density. In addition, as discussed by 
CLMP, eq.(|) can be generalized to multiple streaming as 
a Chapman-Kolmogorov-type relation 

5 h (x,z\M,z f ) = Jdq [l + S^(q\M,z f )] fe[x-x(q, z)] -1, 

with 8d the Dirac function: each fluid element of 
Lagrangian size R Q carries a 'halo density charge' 
rih(M 7 Zf) [1 + 8%(q\M, Zf)] along its trajectory. 

At every point x = q + S, with S the displacement of 
the q-th Lagrangian element smoothed on R a , we assign a 
halo density Sh on the selected mass scale M according to 
eqs.(||) and We extensively use this formulation of our 
bias scheme in § 3 where we test the model locally against 
a high-resolution A-body simulation. 

2.2. The distribution of halos 

In this section we compute the probability distribution 
p(5h) deriving from our bias model. We consider a mildly 
non-linear density field in the laminar regime, though for 
comparisons with simulations we will also adopt the multi- 
strearngeneralization in eq.(^). 

Eq.(fi]) can be recast, using mass conservation, in terms 
of the Jacobian determinant J = ||9x/<9q|| of the map- 
ping from Lagrangian to Eulerian space, q — > x(q, z) = 
q + S(q, z) , namely 1 + 6[x, z] — J(q, z)^ 1 . In the 
Zel'dovich approximation, S(q, z) = — D(z)V q ^ (q), 
where ip is the linear peculiar gravitational potential, 
such that Vq(^o(q) = £o(q). The probability distribu- 
tion of the eigenvalues A Q (q) (a — 1,2,3) of the de- 
formation tensor d 2 <p Q (q)/dq a dqp (Doroshkevich 1970) 
can be used to compute the one-point statistical prop- 
erties of 5h at any redshift z < zt in Eulerian space. 
Let us introduce the variables C a = y5 (fif — 3/^) , and 
Vo = (Ai - 2A 2 + A 3 )/(2 > /(a»? - 3At 2 ) ) , (Reisenegger & 
Miralda-Escude 1995) with /x a (q) ( a — 1,2,3) the invari- 
ants of the deformation tensor. Unlike the original eigen- 
values, these variables are independent, 
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Fig. 1. — The one-point distribution function of 5^. The short 
dashed lines refer to the data obtained by performing 10 5 random 
generations of the variables e ,C and V a (ko = 0.180 h Mpc~ 1 , 
(Tq = 0.217); the filled dots show the TV-body output. The long 
dashed lines represent the linear bias prediction while the solid lines 
are the outcome of the numerical version of our model described in 
§3.2. Top panels: class S (kf = 0.314 h Mpc" 1 ). Bottom panels: 
class L (k f = 0.360 /iMpc" 1 ). 

where < C < oo and -1/2 < V < 1/2. The Ja- 
cobian now reads J = 1 - D{675e Q + 45L>(£ 2 - 5e 2 ) + 
D 2 [2V5£ 3 Vo(3 - 4V 2 ) - 15£ 2 e + 25e 3 ]}/675; it enters in 
eq.(l) through mass conservation. 

The probability p{Sh) can then be computed by Monte 
Carlo generating realizations from the distribution in 
eq.(||). Since we are interested in the Eulerian probability 
and eq.(|q) gives a Lagrangian distribution, we compute 
P($h) =Jd5 p L (S h ,S)/(l + S), where p L (S h ,5) is the joint 
Lagrangian probability for the Eulerian halo and mass 
overdensity fields [cf. eq.(14) in Kofman et al. 1994]. In 
practice, p{o"h) is obtained by: i) generating realizations 
for eo,Co,V ] ii) computing J, and Sh through eq.(0) or 
eq.(|l|); Hi) weighting the contribution to the probability of 
Sh by the factor J. 

3. TESTING THE MODEL 

3.1. Comparing the probability distribution function 

We test our predictions for p(Sh) against a high- 
resolution A-body simulation from the data bank of cos- 
mological simulations provided by the Hydra Consortium 
and produced using the Hydra A-body code (Couchman, 
Thomas & Pierce 1995). The simulation (RUN 501) 
evolves 128 3 particles on a 128 3 cubic mesh with periodic 
boundaries. The box size is 100/i _1 Mpc and the parti- 
cle mass 1.32 x 10 11 /i _1 M . The initial conditions are 
Gaussian with a Cold Dark Matter spectrum with shape 
parameter P = 0.25, density parameter fi = 1 and zero 
cosmological constant. The simulation output corresponds 
to as = 0.64, where as is the rms linear mass fluctua- 
tion in spheres of 8/i. _1 Mpc. At this epoch, the charac- 
teristic virializing halo mass, M», defined by a M — S c , 
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is 0.66 x lO 13 ft- _1 M , i.e. 50.13 particles. To compare 
our predictions to the TV-body outcome we need a halo 
catalog from the simulation. We adopt the spherical over- 
density (SO) halo-finder (Lacey & Cole 1994) to identify 
spherical regions with mean overdensity k = 178, leading 
to 5025 halos with more than 20 particles. We then con- 
sider two classes of objects: class S contains halos with 
0.5 < M/M*< 0.7; class L has 3 < M/M«< 6 (see Ta- 
ble 1). In Figure 1 we plot p(8h) obtained with our bias 
scheme against the iV-body outcomes. The model pre- 
diction for z = Zf = is computed as in §2. a a has 
been tuned to optimize the agreement with the numeri- 
cal outputs. The simulation probability distribution has 
been extracted after smoothing the halo distribution by 
a Gaussian filter W = exp(— k 2 /2k"j) with resolution kf. 
The prediction of a linear Eulerian bias model is also plot- 
ted in Figure 1: the mass distribution in the simulation 
has been smoothed with the same filter used for the halo 
overdensity and the resulting 5 is multiplied by the Eule- 
rian bias factor &mw = 1 + [<Wo^ — 1/^c] (MW), reported 
in Table 1. Our model accurately reproduces the tail of 
the distribution for positive Sh, while for 5h < it favours 
moderate undcrdensitics (5h ~ —0.5) with respect to the 
TV-body simulation. The linear bias prescription instead 
produces much a less skewed distribution with a higher 
peak and severely underestimates the probability of very 
underdense regions. Our model can be further improved 
by adopting the multi-stream version introduced in § 3.2, 
whose predictions, also plotted in Figure 1, are in excellent 
agreement with the TV-body outputs. 

3.2. Cross- correlations 

We also performed a much more severe point-by-point 
test, implementing a fully numerical version of our bias 
scheme as follows. We consider a computational box as 
large as the TV-body one, but sampled with lower resolu- 
tion: 64 3 particles on a 64 3 grid (using 128 3 particles on a 
128 3 mesh gave identical results). Each particle is moved 
to its final position according to the 'truncated' Zel'dovich 
approximation (Coles et al. 1993), that is, prior comput- 
ing the displacement, we remove initial power in high fre- 
quency modes by a Gaussian filter with resolution k a ; we 
require the amplitudes and phases of the linear density 
field to be identical to the simulation ones (at least for the 
Fourier modes present in both grids). Each particle is then 
associated to the linear field e (q) = — V q • S(q). In such 
a way, every particle is endowed with its own halo-density 
charge rih [l + <5^ (q)], computed as described in §2.1. Note 
that, while for a 2 , we used Gaussian smoothing, a M is cal- 
culated with top-hat filtering. The halo-density charge is 
carried by the particles and eventually assigned to the 64 3 
grid through the Triangular Shaped Cloud scheme. The 
corresponding halo overdensity field in the multi-stream 
regime <5™ od is then computed and smoothed with a Gaus- 

Table 1 
Halo parameters 



Class 


TV min 


■^V max 


#halos 
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25 


35 


1021 


0.90 


L 


150 


300 
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Fig. 2. — Cross correlation coefficients vs. filtering scale. Filled 
and empty dots represent, respectively, our model and any Eulerian 
linear bias scheme. 



sian filter with kf < k . The result has to be com- 
pared with the halo field <5| lm extracted from the N-body 
simulation, smoothed on |he same scale. To quantify 
the agreement between the two, we compute their cross- 
correlation coefficient X = (8h° d Sf m )/((J^ od cr^ m ) where 
a % h = (Sh 2 ) 1 / 2 , and the average is performed over the grid 
points. A value |x| = 1 means that the two fields are pro- 
portional, while x — for uncorrelated fields. We tune the 
truncation k Q to optimize our bias scheme. Since x turns 
out to depend very weakly on k (in a wide range around 
its optimal value), we can choose it so that cr™ od = cr" m 
while keeping the maximum allowed value of X - F° r eacn 
resolution this is obtained with a 2 = 1.65 a 2 + C, with 
C = 1.45 for class S and C = 0.54 for class L. Note that 
in the MW model the value of er™ od /V^ lm is not adjustable 
and generally differs from 1. For the resolutions considered 
in Figure 3, o-™ od /of m is 0.81 for class S and 0.90 for class 
L. Moreover, considering smaller smoothing lenghts, this 
ratio can differ from 1 even by 30%. This implies the MW 
method is unable to accurately predict the average prop- 
erties of the bias distribution. Figure 2 reports \ values 
for the optimized model vs. the resolution kf. 

To evaluate the performance of any Eulerian linear bias 
model, we also show the cross-correlation coefficient when 
the smoothed mass density field from the simulation is 
used instead of the output mass of our model. This test 
does not depend on the value of the linear bias. For each 
smoothing length our model reproduces the halo density 
field of the simulation much more accurately than the lin- 
ear biasing scheme. The performance of the two models is 
similar only for very large smoothing lengths. 

In Figure 3 we show the scatter obtained by plotting 
jmod ys £sim j? or our moc iel an d f or the linear bias scheme. 
In this case we adopt the MW bias factor in Table 1. Even 
though on average our bias scheme gives better predictions 
(especially for class S and for underdense regions), some 
scatter in the relation <5™ od vs. <5^ lm persists. To test whe- 
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Fig. 3. — Model predictions for 5^ vs. ./V-body outcome. Left 
panels: our bias scheme (Conserved Halo-Density Charge, CHC). 
Right panels: Eulerian Linear Bias model (ELB) in the MW ver- 
sion. Top panels: class 5 (with kf = 0.260 hMpc -1 ). Bottom pan- 
els: class L (with kf = 0.301 /iMpc -1 ). A l-in-8 random selection 
is shown. 

ther this is due to our simplified dynamics, we generate 
new halo maps taking the particle displacements directly 
from the TV-body simulation. The results are in excel- 
lent agreement with those obtained with the Zel'dovich 
approximation, indicating that the local PS approach, be- 
ing unable to accurately model the Lagrangian halo count- 
ing (Fig. 8 in White 1996), is actually responsible for the 
scatter. We will address this point in a future work. 

4. DISCUSSION AND CONCLUSIONS 

We devised a simple and fast semi-analytical technique 
which allows to study the spatial distribution of dark mat- 
ter halos in terms of their local number density contrast. 
Our method, which is based on a Lagrangian halo identifi- 



cation algorithm plus the Zel'dovich approximation for the 
matter dynamics, was successfully tested against the dis- 
tribution of halos extracted from a high resolution N-body 
simulation. Possible improvements should go in the direc- 
tion of refining the halo selection criterion in Lagrangian 
space, e.g. using the ellipsoidal collapse model or the peak 
theory as in (Bond & Myers 1996; see also Monaco 1998). 

As stressed by CLMP, our model can be applied to study 
the evolution of galaxy biasing, once the relation between 
the galaxies and the hosting dark matter halos is speci- 
fied (e.g. Matarrese et al. 1997). In particular, defining 
the bias field such that 5h = b 5, from eq.(7) one obtains 
b[x(q, z), z) = l + 5%(q\z f )/[l - J(q, z)}. Tegmark & Pee- 
bles (1998) have recently stressed the importance of the 
asymptotic trend of the bias factor. We can analyze this 
issue in the present context by considering a galaxy popu- 
lation conserved in number after an initial merging phase 
(i.e. for varying z at fixed Zf). In the Einstein-de Sitter 
case we recover the 'debiasing' predicted by linear theory: 
b — > 1 as z - * — 1. Differently, if O < 1, b tends to a 
space-dependent value which generally differs from 1; lin- 
ear theory would predict b — > 1 + (bo — l)/D-i(fl ), as 
z ->■ -1, with b Q = b(z = 0) and £>_i = D(z = -1). 

Our method can be used to analyze the coarse-grained 
statistical properties of galaxies and clusters at various 
rcdshifts, e.g. applying semi-analytical techniques to re- 
late the dark matter halo distribution to that of luminous 
objects like galaxies (e.g. Kauffmann, Nusscr & Steinmetz 
1997). It can be further implemented to generate very 
large mock maps of these objects in our past light-cone, 
a problem made of compelling relevance by the ongoing 
wide-field redshift galaxy surveys like the 2 Degree Field 
Survey (2dF) and the Sloan Digital Sky Survey (SDSS). 
Maps of the X-ray cluster distribution may also be pro- 
duced by the present method. 
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by the Italian MURST. CP thanks support from NASA 
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